\chapter{从最佳逼近入手}
\setcounter{equation}{0}
\setcounter{theorem}{0}
\setcounter{definition}{0}

有限元方法是一种经典的微分方程数值求解理论和方法. 目前一般认为, 在上世
纪 50 年代左右, 工程领域已经开始运用类似的方法求解力学问题. 而在数学领
域, 有限元可以看作是古典变分方法的一个发展和补充. 在我国国内, 冯康先生
于 1965 年发表在《应用数学与计算数学》第 2 卷， 第 4 期上的总结性论文
《基于变分原理的差分格式》(全文见 ref 目录下文件 fengkang1965.pdf), 则
普遍被西方学术界认为是中国独立发展有限元方法的一个标志.

作为一种微分方程数值解法, 有限元的最终目的是给出方程的数值解. 为了明确
我们的目标, 突出有限元的基本思想和实现技术, 我们直接跳过全部求解过程,
假设我们已经有了问题的解, 我们应该如何表达它呢? 甚至, 我们可以做的更过
份一点, 假设我们已知问题的解是一个连续函数
\begin{equation}
  f(x), x \in \mathbb{R}^d,
  \label{eq::real_sol_dd}
\end{equation}
这里 $d$ 代表空间的维数. 通常一个物理问题的解, 是在三维空间加上一维时间
中的. 这里我们先考虑其最简形式, $d = 1$, 并不考虑时间. 

既然是数值解, 我们只能用某种逼近的手段表达其函数值的近似值, 或者更准确
地, 我们需要首先确定我们的表达空间, 比如分片线性空间. 对任意的 $f(x)$,
我们需要在表达空间中提供它的最佳逼近元. 当然这里涉及到最佳逼近元的存在
性和构造方法等一系列问题, 而有限元的理论和方法本质上就是对这一系列问题
的一个回答. 我们这里既然跳过了全部中间过程, 直接看结果, 那么先假定我们
目标方程的解析解是一维连续函数 $f(x)$, 而在我们的目标空间, 分片线性空
间中一定存在一个合理的表达 $f_h(x)$.

\section{分片线性空间中的数值解}

设 $\{x_0, x_1, \cdots, x_N\}$ 是区间 $[a, b]$ 上的一个网格剖分. 即
\begin{equation}
  a = x_0 < x_1 < \cdots < x_N = b,
  \label{eq::grid1d}
\end{equation}
构建基函数
\begin{equation}
  \phi_i(x) = \left\{
    \begin{array}{ll}
      \frac{x - x_{i - 1}}{x_i - x_{i - 1}} ,& x \in [x_{i - 1}, x_i], i = 1, 2, \cdots, N; \\
      \frac{x - x_{i + 1}}{x_i - x_{i + 1}} ,& x \in [x_i, x_{i + 1}], i = 0, 1, \cdots, N - 1; \\
      0, & \mbox{其他}.
    \end{array}  
  \right.
  \label{eq::basis_p_lin_1d}
\end{equation}

图 \ref{fig::basis_p_pin_1d} 是基函数 \ref{eq::basis_p_lin_1d} 的示意
图. 这些基函数张成的线性空间即为一维分片线性空间, 记作 $V_h^1$. 也即
\begin{equation}
  V_h^1 = \mbox{span} \{\phi_i, i = 0, 1, \cdots, N\}.
  \label{eq::p_lin_space_1d}
\end{equation}
注意这里称一维是因为 $V_h^1$ 是一维连续函数的分片线性逼近元所构成的线性
空间, 实际上它自身是一个$N + 1$ 维线性空间.

\begin{figure}[!ht]
  \centering
    \begin{tikzpicture}
      \draw[line width=1pt] (0,0) -- (1,0) ;
      \draw[line width=1pt, dashed] (1,0) -- (2,0) ;
      \draw[line width=1pt] (2,0) -- (3,0) -- (4,0);
      \draw[line width=1pt, dashed] (4,0) -- (5,0) ;
      \draw[line width=1pt] (5,0) -- (6,0);
      \filldraw[black] (0,0) circle(2pt) node[anchor=north] {$x_0$} ;
      \filldraw[black] (1,0) circle(2pt) node[anchor=north] {$x_1$} ;
      \filldraw[black] (2,0) circle(2pt) node[anchor=north] {$x_{i - 1}$} ;
      \filldraw[black] (3,0) circle(2pt) node[anchor=north] {$x_{i}$} ;
      \filldraw[black] (4,0) circle(2pt) node[anchor=north] {$x_{i + 1}$} ;
      \filldraw[black] (5,0) circle(2pt) node[anchor=north] {$x_{N - 1}$} ;
      \filldraw[black] (6,0) circle(2pt) node[anchor=north] {$x_{N}$} ;
      \draw[line width=0.5pt] (0, 2) -- (1, 0) ;
      \draw[line width=0.5pt, dashed] (0,0) -- (0, 2) node[anchor=south] {$\phi_0$};    
      \draw[line width=0.5pt] (3, 2) -- (2, 0) ;
      \draw[line width=0.5pt, dashed] (3,0) -- (3, 2) node[anchor=south] {$\phi_i$};    
      \draw[line width=0.5pt] (3, 2) -- (4, 0) ;
      \draw[line width=0.5pt, dashed] (6,0) -- (6, 2) node[anchor=south] {$\phi_N$};    
      \draw[line width=0.5pt] (6, 2) -- (5, 0) ;
    \end{tikzpicture}
    \caption{一维分片线性空间的基函数.}
    \label{fig::basis_p_pin_1d}
\end{figure}

对于 $[a, b]$ 上的一维连续函数 $f(x)$, 在 $V_h^1$ 上, 我们考虑如何构建
它的最佳逼近元
\begin{equation}
  f_h = \sum_{i = 0}^N f_i \phi_i(x),
  \label{eq::p_lin_func_1d}
\end{equation}
也就是说, 我们需要确定 $f_h(x)$ 在网格结点 $x_i$ 上的取值
\begin{equation}
  f_i := f_h(x_i), i = 0, 1, \cdots, N,
  \label{eq::func_val_1d}
\end{equation}
使得 $\|f(x) - f_h(x)\|$ 最小. 这里的范数先取大家熟悉的 L2 范数, 因此
本质上这是一个最佳平方逼近问题, 以后我们会讨论用哪一种范数更加合适. 我
们在数值分析课程中曾经讨论过, $f(x)$ 在 $V_h^1$ 上的最佳平方逼近元
$f_h$ 是它的正交投影, 也即满足
\begin{equation}
  <f_h(x) - f(x), \phi_j> = 0, j = 0, 1, \cdots, N.
  \label{eq::best_ele_L2}
\end{equation}

对于一般的线性空间, 求解该正规方程组是困难的, 除非 $\{\phi_i, i = 0,
1, \cdots, N\}$ 是一组正交基. 但是这里尽管 $\{\phi_i, i = 0,
1, \cdots, N\}$ 并不是两两正交, 但却有局部性,
\begin{equation}
  <\phi_i, \phi_j> = 0, |i - j| \geq 2.
  \label{eq::local_bas}
\end{equation}
或者说, 在 $V_h^1$ 空间, 正规方程组 (\ref{eq::best_ele_L2}) 是三对角方程
组, 其求解是容易的. 

我们用 $K$ 表示正规方程组的系数矩阵，则当 $|i - j| < 2$ 时, 
\begin{equation}
\begin{array}{rcl}
  K_{ij} = <\phi_i, \phi_j> &=& \int_a^b \phi_i(x) \phi_j(x) dx \\
  &=& \sum_{i = 1}^N \int_{x_{i - 1}}^{x_i}\phi_i(x) \phi_j(x) dx \\
  &=& \int_{x_{j - 1}}^{x_j}\phi_i(x) \phi_j(x) dx
      + \int_{x_{j}}^{x_{j + 1}}\phi_i(x) \phi_j(x) dx.    
\end{array}
\label{eq::norm_eq_coeff}
\end{equation}
由积分可加性, 每个子区间 $[x_{i - 1}, x_i]$, 上面通过积分累加给系数矩阵 $K$ 的值为
\begin{equation}
  \begin{array}{rcl}
  c_{00} &=& \int_{x_{i - 1}}^{x_i}\phi_{i - 1}(x) \phi_{i - 1}(x) dx, \\
  c_{01} &=& c_{10} = \int_{x_{i - 1}}^{x_i}\phi_{i - 1}(x) \phi_i(x) dx, \\ 
  c_{11} &=& \int_{x_{i - 1}}^{x_i}\phi_i(x) \phi_i(x) dx.
 \end{array}
  \label{eq::cont_coeff}
\end{equation}
这里 $c_{00}$, $c_{01}$, $c_{10}$ 和 $c_{11}$ 分别累加给
$K_{i - 1, i - 1}$, $K_{i - 1, i}$, $K_{i, i -
  1}$ 和 $K_{ii}$. 一般我们称
\begin{equation}
  C = \left[
    \begin{array}{cc}
      c_{00} & c_{01} \\
      c_{10} & c_{11}
    \end{array}
  \right]
  \label{eq::cont_mat}
\end{equation}
为{\heiti 局部贡献矩阵}. 通过贡献矩阵, 我们可以只遍历一次所有单元 (子区间), 就得
到正规方程组 (\ref{eq::best_ele_L2}) 的系数矩阵 $K$. 而该方程组的右端向
量 $r$ 可以类似地通过积分
$$
r_{i - 1} += \int_{x_{i - 1}}^{x_{i}} f(x) \phi_{i - 1}(x) dx
$$
和
$$
r_i += \int_{x_{i - 1}}^{x_{i}} f(x) \phi_i(x) dx
$$
得到 $e_i$ 单元对它的贡献.

在一维情况, 为计算和编程方便, 我们总是将 $[x_{i - 1}, x_i]$ 变换
到 $[-1, 1]$ 区间, 则局限在其上的基函数为
\begin{equation}
  \psi_0(x) = \frac{1 - x}{2}, \psi_1(x) = \frac{x + 1}{2}.
\end{equation}
这种局限在局部单元上的半个基函数, 通常称为{\heiti 形函数}. 它是基函数
的另一种表示, 只是和单元挂钩, 便于我们局部计算和拼装. 这种变换到一个标
准单元的做法, 在有限元计算中是常见的. 这种附加了形函数信息的标准单元,
被成为{\heiti 模板单元}或{\heiti 参考单元}. 在模板单元上, 除了形函数的
个数事实上是和我们希望的逼近精度有关的, 也就是说 $V_h^1$ 中每个单元上
都是一个分段一次函数, 如果我们要每个单元都是一个分段 $k$ 次多项式, 那
么对应的就是 $V_h^k$ 空间, 该空间的模板单元上, 就会有 $k$ 个形函数, 它
们分布在某种 $k$ 次插值多项式的插值节点上, 这些插值节点, 我们称为
{\heiti 自由度}, 因为本质上我们需要确定的, 就是 $f_h$ 在这些点上的取值,
也就是说是我们问题的实际未知量. 区间 $[a, b]$ 上全部 (由模板单元映射回
的) 自由度形成了正规方程组的完整的未知向量 $f_h$. 有了模板单元, 我们全
部的积分就可以在模板单元上通过数值积分完成.

具体而言, 首先构建模板单元 $e = [-1, 1]$, 自由度共两个, 分别为
$$
\xi_0 = -1, \xi_1 = 1.
$$
对应的形状函数为
$$
\psi_0 = \frac{1 - x}{2}, \psi_1 = \frac{x + 1}{2}.
$$

对一个具体的全局单元 $e_i = [x_{i - 1}, x_i]$, $i = 1, 2, \cdots, N$,
它到模板单元 $[\xi_0, \xi_1]$ 的线性坐标变换为:
$$
[x_{i - 1}, x_i] \mapsto [\xi_0, \xi_1], \xi = (x - x_{i - 1})
   \frac{\xi_1 - \xi_0}{x_{i - 1} - x_i} + \xi_0.
$$
因此对其上的积分有
$$
   \int_{x_{i - 1}}^{x_i} g(x) dx = \int_{x(\xi_0)}^{x(\xi_1)} g(x(\xi)) dx(\xi)
   = \frac{x_{i - 1} - x_i}{\xi_1 - \xi_0} \int_{\xi_0}^{\xi_1} g(x(\xi)) d \xi. 
$$
这里 $g(x)$ 是一般函数, 也即只要在标准单元上计算 $g(x)$ 在 $x
\mapsto \xi$ 这个线性变换下的积分, 然后再乘以变换的比例系数即得原单元上的积分 (贡献).

在数值积分时,
$$
   \int_{\xi_0}^{\xi_1} g(x(\xi)) d\xi \approx \sum_{l = 0}^k g(x(\xi_l)) \omega_{l},
$$
   
其中 $\xi_l$ 和 $\omega_l$ 为对应的高斯积分点和积分权重. 也就是说我们
只需要在模板单元上完成数值积分, 然后将积分值乘以坐标变换比例系数即可.
为此, 我们只需要保存模板单元上的各阶积分高斯积分信息.

注意, 在手算时, 如果做精确积分, 不必做坐标变换, 直接在各单元积分即可.
模板单元只是为了编程方便.

这里对系数矩阵 $K$ 的贡献的形式为 $<\phi_i, \phi_j>$, 因为全部基函数都
是多项式, 因此只要提供足够代数精度的数值积分公式就能得到精确的结果. 而
对于右端项, 我们不能期待 $f$ 总是多项式, 因此这里会产生积分误差, 而对
于逼近问题而言, 这种误差是本质的,不可避免的, 而且已经是最小的.

而更一般的, 如果我们希望数值解有更高的精度, 那么我们可以在 $V_h^k$ 空
间, 也就是分片 $k$ 次多项式空间中构建 $f(x)$ 的最佳逼近元. 显然, 对于
$V_h^1$ 空间中的数值解, 在每个子区间内部, 只有一次精度, 而区间之间, 只
有函数连续性, 导数是间断的. 现在我们将一个子区间 $[x_{i - 1}, x_i]$ 记
作一个单元 $e_i$, $i = 1, 2, \cdots, N$. 那么对于 $V_h^k$ 空间中的数值
解, 我们的任务是在每一个单元内部, 构建一个目标函数 $f$ 的最佳逼近元的
$k$ 次多项式插值, 而在单元之间, 我们暂时仍只要求连续. 比如以 $k = 2$
为例, 为了在单元 $e_i = [x_{i - 1}, x_i]$ 上构建逼近连续目标函数
$f(x)$ 的 $2$ 次多项式表示, 我们可以参考多项式插值节点:
 \begin{equation}
   d_i^{(0)} = x_{i - 1}, d_i^{(1)}
   = \frac{x_{i - 1} + x_i }{2}, d_i^{(2)} = x_i.
   \label{eq::dof_2_1d}
 \end{equation}
 这是一种常见而自然的多项式表示方法. 而在这个三个插值节点上, 我们容易构建 $2$ 次插
 值基函数:
 \begin{equation}
   \begin{array}{rcl}
     \phi_i^{(0)} &=& \frac{(x - d_i^{(1)})(x - d_i^{(2)})}
                      {(d_i^{(0)} - d_i^{(1)})(d_i^{(0)} - d_i^{(2)})}; \\
     \phi_i^{(1)} &=& \frac{(x - d_i^{(0)})(x - d_i^{(2)})}
                      {(d_i^{(1)} - d_i^{(2)})(d_i^{(1)} - d_i^{(2)})}; \\
     \phi_i^{(2)} &=& \frac{(x - d_i^{(0)})(x - d_i^{(1)})}
                      {(d_i^{(2)} - d_i^{(0)})(d_i^{(2)} - d_i^{(1)})}.
   \end{array}
   \label{eq::basis_2_1d}
 \end{equation}
 而此时, 整体逼近问题的对应正规方程组系数矩阵 $K$, 仍然是局部的, 它即
 便不是带状对角矩阵, 但确实是稀疏的. 容易通过数值手段求解.
 
 综合起来, 在有限元框架下, 对一个已知函数的逼近, 我们采用的数值手段是
 先构建 $V_h^k$ 空间, 这里首先确定的是单元剖分. 这种剖分可以不必是等距
 的. 因为本质上我们只需要在每个单元上积分. 然后是单元基函数, 其实我们
 只需要确定模板单元, 然后确定模板单元上的局部形函数即可. 和局部形函数
 一一对应的是自由度的坐标. 最后为了能够有正确的计算精度, 我们必须要提
 供合理的积分公式, 也就是积分的代数精度至少要保证基函数之间的运算结果
 是精确. 这里我们一般会采用高斯积分公式, 可以通过外部手段预先将模板单
 元的高斯积分节点和权重存储待查. 有了这些准备, 我们就可以拼装出正规方
 程组并求得目标函数的最佳逼近. 数值结果以 $f_h$ 在自由度节点上的取值的
 方式存储, 在绘图等实际工作中, 通过插值多项式得到任意点的函数值.

\begin{example} 
  手算 $\sin(x)$, $x \in [0, \pi]$ 在 $V_h^1$ 中的最佳逼近元, 其中
  $$
  x_0 = 0, x_1 = \frac{\pi}{4}, x_2 = \frac{\pi}{3},
  x_3 = \frac{\pi}{2}, x_4 = \frac{2\pi}{3}, x_5 = \pi.
  $$
  
  由于是手算, 我们直接写出全部基函数:
  $$
  \begin{array}{rcl}
    \phi_0(x) &=& 1 - \frac{4}{\pi} x, x \in [0, \frac{\pi}{4}] \\
    \phi_1(x) &=&
    \left \{
    \begin{array}{ll}
      \frac{4}{\pi} x ,& x \in [0, \frac{\pi}{4}]\\
      4 - \frac{12}{\pi} x,& x \in [\frac{\pi}{4}, \frac{\pi}{3}]
    \end{array}
    \right. \\
    \phi_2(x) &=&
    \left \{
    \begin{array}{ll}
      \frac{12}{\pi} x - 3,& x \in [\frac{\pi}{4}, \frac{\pi}{3}]\\
      3 - \frac{6}{\pi} x,& x \in [\frac{\pi}{3}, \frac{\pi}{2}]
    \end{array}
    \right. \\
    \phi_3(x) &=&
    \left \{
    \begin{array}{ll}
      \frac{6}{\pi} x - 2,& x \in [\frac{\pi}{3}, \frac{\pi}{2}]\\
      4 - \frac{6}{\pi} x,& x \in [\frac{\pi}{2}, \frac{2 \pi}{3}]
    \end{array}
    \right. \\
    \phi_4(x) &=&
    \left \{
    \begin{array}{ll}
      \frac{6}{\pi} x - 3,& x \in [\frac{\pi}{2}, \frac{2 \pi}{3}]\\
      3 - \frac{3}{\pi} x,& x \in [\frac{2 \pi}{3}, \pi]
    \end{array}
    \right. \\
    \phi_5(x) &=& \frac{3}{\pi} x - 2, x \in [\frac{2 \pi}{3}, \pi]. 
  \end{array}
  $$
  以上基函数在定于域外全为零.

  然后就能直接通过积分 (\ref{eq::norm_eq_coeff}) 得到正规方程组系数矩
  阵 $K$ 和右端向量 $r$. 为了方便计算,考虑到以上基函数都是分段的, 我们
  利用积分可加性, 按单元次序逐个拼装矩阵 $K$ 和 $r$.

  对 $e_1 = [x_0, x_1] = [0, \frac{\pi}{4}]$, 其上积分不为零的基函数为
  $\phi_0(x)$ 和 $\phi_1(x)$, 对系数矩阵 $K$ 的局部贡献为:
  $$
  \begin{array}{rcl}
    C^{(1)}_{00} &=& \int_0^{\frac{\pi}{4}} \phi_0^2(x) dx
    = \int_0^{\frac{\pi}{4}} \left(1 - \frac{4}{\pi} x\right)^2 dx
    = \frac{\pi}{12} \approx 0.2618, \\
    C^{(1)}_{01} = C^{(1)}_{10} &=& \int_0^{\frac{\pi}{4}} \phi_0(x) \phi_1(x) dx
    = \int_0^{\frac{\pi}{4}} \left(1 - \frac{4}{\pi} x\right) \frac{4}{\pi} x dx
    = \frac{\pi}{24} \approx 0.1309, \\
    C^{(1)}_{11} &=& \int_0^{\frac{\pi}{4}} \phi_1^2(x) dx
    = \int_0^{\frac{\pi}{4}} \left(\frac{4}{\pi} x\right)^2 dx
    = \frac{\pi}{12} \approx 0.2618.
  \end{array}
  $$
  而对右端项 $r$ 的局部贡献为:
  $$
  \begin{array}{rcl}
    c^{(1)}_0 &=& \int_0^{\frac{\pi}{4}} \sin(x) \phi_0(x) dx
    = \int_0^{\frac{\pi}{4}} \sin(x) \left(1 - \frac{4}{\pi} x\right) dx
    \approx 0.09968 , \\
    c^{(1)}_1 &=& \int_0^{\frac{\pi}{4}} \sin(x) \phi_1(x) dx
    = \int_0^{\frac{\pi}{4}} \sin(x) \frac{4}{\pi} x dx
    \approx 0.1932.
  \end{array}
  $$
  
  这里右上标 $^{(1)}$ 表示是 $e_1$ 单元的贡献. 类似地, $e_2$, $e_3$,
  $e_4$ 和 $e_5$ 对应的局部贡献为:
  $$
  \begin{array}{llll}
    e_2:& C_{11}^{(2)} \approx 0.08727,& C_{12}^{(2)} = C_{21}^{(2)} \approx 0.04363,
    & C_{22}^{(2)} \approx 0.08727,\\
    & c_1^{(2)} \approx 0.1001,& c_2^{(2)} \approx 0.1070,&\\ 
    e_3:& C_{22}^{(3)} \approx 0.1745,& C_{23}^{(3)} = C_{32}^{(3)} \approx 0.08727,
    & C_{33}^{(3)} \approx 0.1745,\\
    & c_2^{(3)} \approx 0.2441,& c_3^{(3)} \approx 0.2559,&\\ 
    e_4:& C_{33}^{(4)} \approx 0.1745,& C_{34}^{(4)} = C_{43}^{(4)} \approx 0.08727,
    & C_{44}^{(4)} \approx 0.1745,\\
    & c_3^{(4)} \approx 0.2559,& c_4^{(4)} \approx 0.2441,&\\ 
    e_5:& C_{44}^{(5)} \approx 0.3491,& C_{45}^{(5)} = C_{54}^{(5)} \approx 0.1745,
    & C_{55}^{(5)} \approx 0.3491,\\
    & c_4^{(5)} \approx 0.3270,& c_5^{(5)} \approx 0.1730. & 
  \end{array}
  $$

  而整体系数矩阵 $K$ 和右端向量 $r$ 由上述局部贡献拼装而成:
  $$
  \begin{array}{rcl}
  K &=& \left[
    \begin{array}{llllll}
      C_{00}^{(1)} & C_{01}^{(1)} & 0 & 0 & 0 & 0 \\
      C_{10}^{(1)} & C_{11}^{(1)} + C_{11}^{(2)} & C_{12}^{(2)} & 0 & 0 & 0 \\
      0 & C_{21}^{(2)} & C_{22}^{(2)} + C_{22}^{(3)} & C_{23}^{(3)} & 0 & 0 \\
      0 & 0 & C_{32}^{(3)} & C_{33}^{(3)} + C_{33}^{(4)} & C_{34}^{(4)} & 0 \\
      0 & 0 & 0 & C_{43}^{(4)} & C_{44}^{(4)} + C_{44}^{(5)} & C_{45}^{(5)} \\
      0 & 0 & 0 & 0 & C_{54}^{(5)} & C_{55}^{(5)} 
    \end{array}
  \right]\\
  &\approx& \left[
    \begin{array}{llllll}
      0.2618 & 0.1309 & 0 & 0 & 0 & 0 \\
      0.1309 & 0.3491 & 0.04363 & 0 & 0 & 0 \\
      0 & 0.04363 & 0.2618 & 0.08727 & 0 & 0 \\
      0 & 0 & 0.08727 & 0.3491 & 0.08727 & 0 \\
      0 & 0 & 0 & 0.08727 & 0.5236 & 0.1745 \\
      0 & 0 & 0 & 0 & 0.1745 & 0.3491
    \end{array}
  \right], \\
  r &=& \left[c_{0}^{(1)}, c_{1}^{(1)} + c_{1}^{(2)}, c_{2}^{(2)} + c_2^{(3)},
    c_3^{(3)} + c_3^{(4)}, c_4^{(4)} + c_4^{(5)}, c_5^{(5)}\right]^T \\
  &\approx& \left[0.09968, 0.2933, 0.3511, 0.5118, 0.5711, 0.1730\right]^T.
  \end{array}
  $$

  求解线性方程组
  $$
  K f_h = r,
  $$
  可得 $\sin(x)$ 在 $V_h^1$ 空间的最佳逼近元
  $$
  f = [0.0193, 0.7228, 0.8809, 1.0192, 0.9067, 0.0423]^T.
  $$

  图 \ref{fig::sin_V1} 显示了实际数值结果.
  \label{ex::sin_V1_hand}
\end{example}
   
\begin{figure}[!ht]
  \centering
  \includegraphics[width=0.75\textwidth]{images/sin_V1.eps}
  \caption{$\sin(x)$ 在 $V_h^1$ 上的最佳逼近, 虚线是 $\sin(x)$, 实折线
    是逼近结果. $x$ 轴上的黑点是自由度的坐标, 自由度垂直向上的点虚线是
    对应基函数的取值, 垂直点虚线上方的黑点是 $f_h$ 在每个自由度的取值.
    这是一个分片线性的结果.}
  \label{fig::sin_V1}
\end{figure}

\section{程序实现}

现在来考虑对例 \ref{ex::sin_V1_hand} 这样的最佳逼近问题, 如何通过程序
自动计算.  在本书中, 我们将采用 C++ 实现全部算法.

从之前的计算我们已经明确, 这里最本质的计算分为两个部分:

\begin{enumerate}
\item 拼装系数矩阵 $K$ 和右端向量 $r$;
\item 线性方程组求解.  
\end{enumerate}

在第一步, 我们需要遍历全部单元, 然后通过数值积分计算每个单元对系数矩阵
和右端向量的贡献. 这里为了统一积分格式, 我们需要有一个模板单元. 目前我
们暂时不考虑第二步如何实现, 直接采用一个熟悉的线性方程组求解算法, 如共
轭梯度法. 注意, 由于系数矩阵是稀疏的, 所以这里必须采用某种迭代法来解线
性方程组. 而直接法将导致矩阵稀疏性的丧失 (除非是三对角矩阵这样的特殊情
况).

在明确了算法流程之后, 我们会注意到有限元程序对数据结构有很强的要求. 几
何网格, 模板单元, 自由度, 基函数, 高斯积分信息, 稀疏矩阵, 向量, 数值解,等
等. 这些都需要合理规划. 这里首先要划分结构化数据和非结构化数据. 所谓结
构化数据, 比如常见的{\heiti 结构化网格}, 指的是每一个网格都可以通过计
算临时得到, 而不需存储实际的网格信息. 比如如果是 $[a, b]$ 区间的等距网
格, 那么就有规律
$$
x_i = a + i h, h = \frac{b - a}{N}, i = 0, 1, \cdots, N.
$$
因此不用存储. 但我们在例 \ref{ex::sin_V1_hand} 中, 则不是这样, 因为我们的
网格划分是不等距的, 所以我们必须存储每一个划分节点的实际坐标. 这种就是典
型的{\heiti 无结构网格}.

结构化数据的优点是高效, 计算方便, 节省内存. 但是无结构网格具有更好的几
何匹配性, 更加贴近实际问题, 特别是在 2 维和 3 维.  一个基本原则是, 如
果在不损失本质精度的前提下, 能够应用结构化数据, 那么就应该采用. 但是在
很多实际应用问题中, 无结构数据几乎是必须的. 因此我们两种实现方式, 都应
该掌握.

在 finite-element-method 目录下, 提供的 example1.cpp 是例
\ref{ex::sin_V1_hand} 的一个实现. 它需要数值代数包 Eigen 提供的代数运算支持.
这是一个标准的包, 在各 Linux 的发行版中都可以找到. 编译运行之后, 会得
到 exampl1.m, 在 Matlab 或者 Octave 中运行, 就能得到和图 \ref{fig::sin_V1}
类似的结果.

注意这个 example1.cpp 只是一个练习性质的代码, 它连同它的数据结构类, 都
还远远称不上完美. 分析, 讨论和改进这些代码, 也是我们这门课程的主要内容
之一.


